function [ Theta1,Theta2 ] = nnretrain(X,y,initial_Theta1,initial_Theta2,lambda,maxiter)
%    compute (input->hidden) theta1 and (hidden->output)theta2
%   [ Theta1,Theta2 ] = NNretrain(X,y,initial_Theta1,initial_Theta2,lambda)
%

if nargin==3 || isempty(initial_Theta2)
    initial_Theta2 = initial_Theta1.theta2;
    lambda = initial_Theta1.lambda;
    initial_Theta1 = initial_Theta1.theta1;
end


input_layer_size  = size(X,2);  
hidden_layer_size = size(initial_Theta1,1);
num_labels = size(initial_Theta2,1);


m = size(X, 1);

if nargin<5
    lambda = 0;
end

if nargin<6
    maxiter = 10;
end

initial_nn_params = [initial_Theta1(:) ; initial_Theta2(:)];
 

options = optimset('MaxIter', maxiter);

costFunction = @(p) nnCostFunction(p, ...
                                   input_layer_size, ...
                                   hidden_layer_size, ...
                                   num_labels, X, y, lambda);
                                   
                                   
[nn_params, ~] = fmincg(costFunction, initial_nn_params, options);


Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), ...
                 hidden_layer_size, (input_layer_size + 1));

Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), ...
                 num_labels, (hidden_layer_size + 1));
             
if nargout==1
    Theta1.theta1 = Theta1;
    Theta1.theta2 = Theta2;
    Theta1.lambda = lambda;
end

end